%	Generate Debt policy functions and asset pricing measure for a case with 3 states of LANCIA, RUSSO, and WORRALL (2024)   %

clear all; close all;
tic
warning('off')

cd ..\
cd Stored_file

load Policy_3state.mat

cd ..\
cd Matlab_functions

% generate useful policies
eo1=e1-ey1;
eo2=e2-ey2;
eo3=e3-ey3;

v1aut = (utility(ey1,gamma)+betta*(pi1*utility(eo1,gamma)+(pi2)*utility(eo2,gamma)+(pi3)*utility(eo3,gamma))); % intertemporal utility in autarky in state 1
v2aut = (utility(ey2,gamma)+betta*(pi1*utility(eo1,gamma)+(pi2)*utility(eo2,gamma)+(pi3)*utility(eo3,gamma))); % intertemporal utility in autarky in state 2
v3aut = (utility(ey3,gamma)+betta*(pi1*utility(eo1,gamma)+(pi2)*utility(eo2,gamma)+(pi3)*utility(eo3,gamma))); % intertemporal utility in autarky in state 2
waut = pi1*utility(eo1,gamma)+(pi2)*utility(eo2,gamma)+(pi3)*utility(eo3,gamma);
wmax=max(wgrid);

% restrict on the relevant state
w_thr = max(wgrid(((P_new(1)-P_new<10^(-10)))));
wgrid_new = [w_thr:0.001:max(wgrid)];
cy1=interp1(wgrid,cy1,wgrid_new,'spline');
cy2=interp1(wgrid,cy2,wgrid_new,'spline');
cy3=interp1(wgrid,cy3,wgrid_new,'spline');
w1=interp1(wgrid,w1,wgrid_new,'spline');
w2=interp1(wgrid,w2,wgrid_new,'spline');
w3=interp1(wgrid,w3,wgrid_new,'spline');
wgrid = wgrid_new;

figure('name','promises')
plot(wgrid,w1,'color',[200 200 200]/255,'LineWidth',1)
hold on
plot(wgrid,w2,'color', [128 128 128]/255,'LineWidth',1)
hold on
plot(wgrid,w3,'k','LineWidth',1)
hold on
plot(wgrid,wgrid,'k:','LineWidth',1)
hold on

% next period consumption
cy11=interp1(wgrid,cy1,w1,'spline');
cy12=interp1(wgrid,cy2,w1,'spline');
cy13=interp1(wgrid,cy3,w1,'spline');
cy21=interp1(wgrid,cy1,w2,'spline');
cy22=interp1(wgrid,cy2,w2,'spline');
cy23=interp1(wgrid,cy3,w2,'spline');
cy31=interp1(wgrid,cy1,w3,'spline');
cy32=interp1(wgrid,cy2,w3,'spline');
cy33=interp1(wgrid,cy3,w3,'spline');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%% FIRST BEST ALLOCATION %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
b1fb=(ey1-delta/(delta+betta))./ey1;
b2fb=(ey2-delta/(delta+betta))./ey2;
b3fb=(ey3-delta/(delta+betta))./ey3;

cy1fb=delta/(delta+betta);
cy2fb=delta/(delta+betta);
cy3fb=delta/(delta+betta);

phiy1fb = (cy1fb./ey1);
phiy2fb = (cy2fb./ey2);
phiy3fb = (cy3fb./ey3);

phio1fb = ((1-cy1fb)./ey1);
phio2fb = ((1-cy2fb)./ey2);
phio3fb = ((1-cy3fb)./ey3);

q11fb = pi1.*betta.*(phiy1fb./phio1fb);
q21fb = pi1.*betta.*(phiy2fb./phio1fb);
q31fb = pi1.*betta.*(phiy3fb./phio1fb);
q12fb = pi2.*betta.*(phiy1fb./phio2fb);
q22fb = pi2.*betta.*(phiy2fb./phio2fb);
q32fb = pi2.*betta.*(phiy3fb./phio2fb);
q13fb = pi3.*betta.*(phiy1fb./phio3fb);
q23fb = pi3.*betta.*(phiy2fb./phio3fb);
q33fb = pi3.*betta.*(phiy3fb./phio3fb);

% SDF
m11fbA = delta.*(phiy1fb./phiy1fb);
m21fbA = delta.*(phiy2fb./phiy1fb);
m31fbA = delta.*(phiy3fb./phiy1fb);
m12fbA = delta.*(phiy1fb./phiy2fb);
m22fbA = delta.*(phiy2fb./phiy2fb);
m32fbA = delta.*(phiy3fb./phiy2fb);
m13fbA = delta.*(phiy1fb./phiy3fb);
m23fbA = delta.*(phiy2fb./phiy3fb);
m33fbA = delta.*(phiy3fb./phiy3fb);

m11fbB = (betta/delta).*(phiy1fb./phio1fb);
m21fbB = (betta/delta).*(phiy1fb./phio1fb);
m31fbB = (betta/delta).*(phiy1fb./phio1fb);
m12fbB = (betta/delta).*(phiy2fb./phio2fb);
m22fbB = (betta/delta).*(phiy2fb./phio2fb);
m32fbB = (betta/delta).*(phiy2fb./phio2fb);
m13fbB = (betta/delta).*(phiy3fb./phio3fb);
m23fbB = (betta/delta).*(phiy3fb./phio3fb);
m33fbB = (betta/delta).*(phiy3fb./phio3fb);

m11fb = m11fbA*m11fbB;
m21fb = m21fbA*m21fbB;
m31fb = m31fbA*m31fbB;
m12fb = m12fbA*m12fbB;
m22fb = m22fbA*m22fbB;
m32fb = m32fbA*m32fbB;
m13fb = m13fbA*m13fbB;
m23fb = m23fbA*m23fbB;
m33fb = m33fbA*m33fbB;

% first best primary surplus/deficit
tau1fb = b1fb - (q11fb.*b1fb+q12fb.*b2fb+q13fb.*b3fb);
tau2fb = b2fb - (q21fb.*b1fb+q22fb.*b2fb+q23fb.*b3fb);
tau3fb = b3fb - (q31fb.*b1fb+q32fb.*b2fb+q33fb.*b3fb);

TBRfb=delta*(pi1*b1fb+pi2*b2fb+pi3*b3fb);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%% SECOND BEST ALLOCATION %%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% critical level of debt
bcrit=1-exp(-betta*(w_thr-waut));

% state-contingent outstanding bond
b1 = (ey1-cy1)./ey1;
b2 = (ey2-cy2)./ey2;
b3 = (ey3-cy3)./ey3;

% new issued state-contingent bond
b11 = (ey1-cy11)./ey1;
b12 = (ey2-cy12)./ey2;
b13 = (ey3-cy13)./ey3;
b21 = (ey1-cy21)./ey1;
b22 = (ey2-cy22)./ey2;
b23 = (ey3-cy23)./ey3;
b31 = (ey1-cy31)./ey1;
b32 = (ey2-cy32)./ey2;
b33 = (ey3-cy33)./ey3;

plot(b1,b11,'r',b1,b12,'g',b1,b13,'b')

% inverse of marginal utility
phiy1 = (cy1./ey1);
phiy2 = (cy2./ey2);
phiy3 = (cy3./ey3);

phio11 = ((1-cy11)./ey1);
phio12 = ((1-cy12)./ey2);
phio13 = ((1-cy13)./ey3);
phio21 = ((1-cy21)./ey1);
phio22 = ((1-cy22)./ey2);
phio23 = ((1-cy23)./ey3);
phio31 = ((1-cy31)./ey1);
phio32 = ((1-cy32)./ey2);
phio33 = ((1-cy33)./ey3);

phiy11 = ((cy11)./ey1);
phiy12 = ((cy12)./ey2);
phiy13 = ((cy13)./ey3);
phiy21 = ((cy21)./ey1);
phiy22 = ((cy22)./ey2);
phiy23 = ((cy23)./ey3);
phiy31 = ((cy31)./ey1);
phiy32 = ((cy32)./ey2);
phiy33 = ((cy33)./ey3);

% SDF
m11 = betta.*(phiy1./phio11);
m12 = betta.*(phiy1./phio12);
m13 = betta.*(phiy1./phio13);
m21 = betta.*(phiy2./phio21);
m22 = betta.*(phiy2./phio22);
m23 = betta.*(phiy2./phio23);
m31 = betta.*(phiy3./phio31);
m32 = betta.*(phiy3./phio32);
m33 = betta.*(phiy3./phio33);

figure('name','SDF')
plot(b1,m11,'r',b1,m12,'g',b1,m13,'b')

% state price
q11 = pi1.*m11;
q12 = pi2.*m12;
q13 = pi3.*m13;
q21 = pi1.*m21;
q22 = pi2.*m22;
q23 = pi3.*m23;
q31 = pi1.*m31;
q32 = pi2.*m32;
q33 = pi3.*m33;

figure('name','state price')
plot(b1,q11,'r',b1,q12,'b',b1,q13,'g')
hold on

% bond revenue
BR1=q11.*b11;
BR2=q12.*b12;
BR3=q13.*b13;
TBR=BR1+BR2+BR3;
BRmin=interp1(b1,TBR,bcrit,'spline');

figure('name','Bond revenue')
plot(b1,TBR)

% primary surplus
tau = b1 - BR1-BR3-BR3;

%----------------------------------------------
% Outcome Variables
%----------------------------------------------

cd ..\
cd Stored_file
save('Debt_policy_3state.mat', 'b1', 'b11', 'b12', 'b13', 'bcrit', 'b1fb', 'b2fb', 'b3fb','TBR','tau','tau1fb','tau2fb','tau3fb','BRmin')

